import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from matplotlib import gridspec
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import random
from scipy.spatial.distance import pdist, squareform

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# 固定参数
klist = [0.0022, 0.1369, 0.1754, 0.1616, 0.3022]
khlist = [0.2989, 0.6281, 0.8579, 0.7592]
# bound = np.array([[30, 50], [50, 90], [6, 10]])  # 温度, 湿度, 固含量的边界
bound = np.array([[30 + np.random.uniform(0, 2), 50 + np.random.uniform(-2, 0)], [50 + np.random.uniform(0, 2), 90 + np.random.uniform(-1, 0)], [6 + np.random.uniform(0, 1), 10 + np.random.uniform(-1, 0)]])  # 温度, 湿度, 固含量的边界


# 定义孔面积占比计算函数
def f(T, H, SC, klist, khlist):
    k_1, k_3, k_4, k_5, k_6 = klist
    k_H, S_C, S_S, eta_0 = khlist

    solid_content = SC / 100
    T_k = T + 273.15  # 转换为开尔文温度

    # 常量定义
    C_0 = 24e-3
    m_D_0 = 24e-3
    m_S = 6e-3
    m_C = (m_D_0 + m_S) * solid_content
    ro_D, ro_S, ro_C = 0.948e3, 1.261e3, 1.3e3
    V = m_D_0 / ro_D + m_S / ro_S + m_C / ro_C
    A, B, C = 6.09451, 2725.96, 28.209
    A_0 = -k_1 * (10 ** (A - B / (T_k + C)) * 133.322 * (1 - k_H * H / 100) / ro_D / V)
    n_cf = 6
    k_boltz = 1.380649e-23
    r = 0.3413e-9
    N_A = 6.02214076e23

    # 随时间变化的函数

    def m_D(t):
        return C_0 * np.exp(A_0 * t)

    def m_S_out(t):
        return max(0, m_S - m_D(t) * S_S)

    def m_C_out(t):
        return max(0, m_C - m_D(t) * S_C)

    def phi_C_out(t):
        return m_C_out(t) / ro_C / V

    def eta(temp, t):
        return eta_0 * (1 + 2.5 * phi_C_out(t))

    def D(temp, t):
        return k_boltz * temp / (n_cf * np.pi * r * eta(temp, t))

    def v(temp, t):
        return k_3 * (D(temp, t) ** 0.5)

    def n_density(t):
        return m_S_out(t) / (ro_S * 4 / 3 * np.pi * r ** 3) / V

    def Z(temp, t):
        return (2 ** 0.5) * n_density(t) * np.pi * (r ** 2) * v(temp, t)

    # 确定蒸发结束时间
    tf = 0
    for t in np.arange(0, 500, 0.01):
        if m_D(t) < 0.024 * 0.1:
            tf = t
            break

    # 解微分方程
    dy = lambda m_s, t: (k_4 * Z(T_k, t) + k_5 * v(T_k, t) * (m_S_out(t) - m_s) / V * m_s)
    t_points = np.arange(0.1, tf, 0.01)
    sol = odeint(dy, 0, t_points)

    return k_6 * sol[-1][0] if len(sol) > 0 else 0


# 包装目标函数
def objective_func(individual):
    T, H, SC = individual
    return f(T, H, SC, klist, khlist)


# 改进的遗传算法实现
def Improved_GA(objective_func, bounds, pop_size=50, max_generations=200):
    dim = len(bounds)
    lb, ub = bounds[:, 0], bounds[:, 1]

    # 初始化种群
    population = np.random.uniform(lb, ub, (pop_size, dim))
    fitness = np.array([objective_func(ind) for ind in population])

    # 记录优化过程
    history = {
        'best_fitness': [],
        'best_individual': [],
        'avg_fitness': [],
        'diversity': [],
        'temperature': [],
        'humidity': [],
        'solid_content': [],
        'all_populations': [population.copy()],
        'all_fitness': [fitness.copy()]
    }

    # 记录初始状态
    best_idx = np.argmax(fitness)
    best_individual = population[best_idx].copy()
    best_fitness = fitness[best_idx]

    history['best_fitness'].append(best_fitness)
    history['best_individual'].append(best_individual)
    history['avg_fitness'].append(np.mean(fitness))
    history['temperature'].append(best_individual[0])
    history['humidity'].append(best_individual[1])
    history['solid_content'].append(best_individual[2])

    # 计算初始多样性
    diversity = np.mean(np.std(population, axis=0))
    history['diversity'].append(diversity)

    # 遗传算法参数
    crossover_rate = 0.9
    mutation_rate = 0.2
    stagnation_counter = 0
    last_improvement = 0

    for gen in range(max_generations):
        # 动态调整参数 - 基于多样性和停滞情况
        if diversity < 0.1:
            mutation_rate = min(0.5, mutation_rate * 1.2)
            crossover_rate = max(0.7, crossover_rate * 0.9)
        elif stagnation_counter > 10:
            mutation_rate = min(0.5, mutation_rate * 1.5)
            crossover_rate = max(0.6, crossover_rate * 0.8)
        else:
            mutation_rate = max(0.1, mutation_rate * 0.95)
            crossover_rate = min(0.95, crossover_rate * 1.05)

        # 选择 - 锦标赛选择
        selected_indices = []
        tournament_size = 3
        for _ in range(pop_size):
            contestants = np.random.choice(pop_size, tournament_size)
            winner = contestants[np.argmax(fitness[contestants])]
            selected_indices.append(winner)

        selected_population = population[selected_indices]

        # 交叉 - 模拟二进制交叉 (SBX)
        offspring = np.zeros_like(population)
        for i in range(0, pop_size, 2):
            parent1 = selected_population[i]
            parent2 = selected_population[(i + 1) % pop_size]
            child1 = parent1.copy()
            child2 = parent2.copy()

            if random.random() < crossover_rate:
                # 对每个维度进行交叉
                for j in range(dim):
                    if random.random() <= 0.5:
                        # 确保不会交换相同位置
                        if abs(parent1[j] - parent2[j]) > 1e-10:
                            beta = 1.0
                            u = random.random()
                            if u <= 0.5:
                                beta = (2.0 * u) ** (1.0 / (1.0 + 1.0))
                            else:
                                beta = (1.0 / (2.0 * (1.0 - u))) ** (1.0 / (1.0 + 1.0))

                            # 计算子代
                            child1[j] = 0.5 * ((1 + beta) * parent1[j] + (1 - beta) * parent2[j])
                            child2[j] = 0.5 * ((1 - beta) * parent1[j] + (1 + beta) * parent2[j])

                            # 确保在边界内
                            child1[j] = max(lb[j], min(ub[j], child1[j]))
                            child2[j] = max(lb[j], min(ub[j], child2[j]))

            offspring[i] = child1
            offspring[(i + 1) % pop_size] = child2

        # 变异 - 多项式变异
        for i in range(pop_size):
            for j in range(dim):
                if random.random() < mutation_rate:
                    # 计算变异量
                    delta = 0.0
                    u = random.random()
                    if u < 0.5:
                        delta = (2.0 * u) ** (1.0 / (1.0 + 1.0)) - 1.0
                    else:
                        delta = 1.0 - (2.0 * (1.0 - u)) ** (1.0 / (1.0 + 1.0))

                    # 应用变异
                    offspring[i][j] = offspring[i][j] + delta * (ub[j] - lb[j])
                    offspring[i][j] = max(lb[j], min(ub[j], offspring[i][j]))

        # 移民 - 引入新个体保持多样性
        num_immigrants = int(pop_size * 0.1)
        immigrant_indices = np.random.choice(pop_size, num_immigrants, replace=False)
        for idx in immigrant_indices:
            # 在最优解附近扰动产生新个体
            perturbation = np.random.uniform(-0.2, 0.2, dim) * (ub - lb)
            offspring[idx] = best_individual + perturbation
            offspring[idx] = np.clip(offspring[idx], lb, ub)

        # 评估子代适应度
        offspring_fitness = np.array([objective_func(ind) for ind in offspring])

        # 环境选择 - 精英保留策略
        # 合并父代和子代
        combined_population = np.vstack((population, offspring))
        combined_fitness = np.hstack((fitness, offspring_fitness))

        # 选择适应度最高的个体
        sorted_indices = np.argsort(combined_fitness)[::-1]  # 从大到小排序
        next_population = combined_population[sorted_indices[:pop_size]]
        next_fitness = combined_fitness[sorted_indices[:pop_size]]

        # 更新种群
        population = next_population
        fitness = next_fitness

        # 更新最优解
        current_best_idx = np.argmax(fitness)
        current_best_individual = population[current_best_idx]
        current_best_fitness = fitness[current_best_idx]

        improvement = False
        if current_best_fitness > best_fitness:
            best_individual = current_best_individual.copy()
            best_fitness = current_best_fitness
            stagnation_counter = 0
            last_improvement = gen
            improvement = True
        else:
            stagnation_counter += 1

        # 计算多样性
        diversity = np.mean(np.std(population, axis=0))

        # 记录当前代信息
        history['best_fitness'].append(best_fitness)
        history['best_individual'].append(best_individual.copy())
        history['avg_fitness'].append(np.mean(fitness))
        history['diversity'].append(diversity)
        history['temperature'].append(best_individual[0])
        history['humidity'].append(best_individual[1])
        history['solid_content'].append(best_individual[2])
        history['all_populations'].append(population.copy())
        history['all_fitness'].append(fitness.copy())

        print(
            f'代 {gen + 1}/{max_generations}, 最佳孔面积占比: {best_fitness:.4f}%, 多样性: {diversity:.4f}, 停滞: {stagnation_counter}, 变异率: {mutation_rate:.3f}')

    return best_individual, best_fitness, history


# 运行改进的遗传算法优化
best_individual, best_fitness, history = Improved_GA(objective_func, bound, pop_size=60, max_generations=20)


# 可视化优化过程
def plot_GA_history(history):
    plt.figure(figsize=(16, 12))
    plt.suptitle('多孔膜制备条件优化过程 - 遗传算法', fontsize=16)

    # 子图1：损失函数变化
    plt.subplot(2, 2, 1)
    plt.semilogy(history['best_fitness'], 'r-', label='最佳适应度')
    plt.semilogy(history['avg_fitness'], 'b--', label='平均适应度')
    plt.xlabel('代数')
    plt.ylabel('孔面积占比 (%)')
    plt.title('孔面积占比优化过程')
    plt.legend()
    plt.grid(True, which="both", ls="--")

    # 子图2：温度衰减曲线
    plt.subplot(2, 2, 2)
    plt.plot(history['temperature'], 'g-')
    plt.plot(history['humidity'], 'b-')
    plt.plot(history['solid_content'], 'r-')
    plt.xlabel('代次')
    plt.ylabel('温度、湿度、固含量')
    plt.title('制备条件曲线')
    plt.grid(True, ls="--")

    # 子图2：温度衰减曲线
    plt.subplot(2, 2, 3)
    plt.plot(history['diversity'], 'm-')
    plt.xlabel('代数')
    plt.ylabel('多样性指标')
    plt.title('群体多样性变化')
    plt.grid(True, ls="--")

    plt.tight_layout()
    plt.show()


# 绘制优化过程
plot_GA_history(history)

# 输出最优解
print("\n最优制备条件:")
print(f"温度: {best_individual[0]:.2f} °C")
print(f"湿度: {best_individual[1]:.2f} %")
print(f"固含量: {best_individual[2]:.2f} %")
print(f"最大孔面积占比: {best_fitness:.4f} %")
